function handles = mapROImeans(handles, figFlag)
% mapROImeans
% handles = mapROImeans(handles, figFlag)
%
% old description:
% MAPROIMEANS   Returns row, column, or overall means and errors of a mapStack's ROI.
% 	Syntax: [vecX, vecY, vecAvg, vecErr] = mapROImeans(mapStack, roi, option, figFglag, Xvec)
%   ROI is a 4 element vector: [rowStart rowEnd colStart colEnd]
%       If roi = [], roi is set to full size of matrix.
%   OPTIONS include: 
%       'row' -- average along rows; i.e., 'side' or Y view
%       'col' -- average along columns; i.e., 'bottom' or X view
%       'all' -- average twice to get overall mean within ROI
% 	XVEC is a user supplied vector for plotting; i.e., it is the 'XData'. Use [] as default.
%   The algorithm is to calculate the average over the ROI first for each individual
%       map, then to obtain the overall mean and SEM.
% 
% 	All the maps in the stack will be evaluated (i.e., no ROI in Z).
%   The averaging is performed using nanmean, not mean.
% 
%   Example:
%   [vecX, vecY, vecAvg, vecErr] = mapROImeans(M, [1 16 1 16], 'col', 1, [])
%
%   To analyze multiple ROIs as one, delete the unwanted data from the
%   map stack, then call this function with roi = []. See mapROIplotsNL for
%   examples.
%
%   See also mapAverage, mapROIplotsNL.

% Changes: gs Mar 2004 - basics
% ***********************************************************************************************
if nargin == 1
    figFlag = 1;
end

mapStack = handles.analysis.mapAverager.mapStack;
option = handles.analysis.mapAverager.ROI.option;
switch option
    case 'col'
        rowStart = handles.analysis.mapAverager.ROI.colROIfirstRow;
        rowEnd = handles.analysis.mapAverager.ROI.colROIlastRow;
        colStart = handles.analysis.mapAverager.ROI.colROIfirstCol;
        colEnd = handles.analysis.mapAverager.ROI.colROIlastCol;
    case 'row'
        rowStart = handles.analysis.mapAverager.ROI.rowROIfirstRow;
        rowEnd = handles.analysis.mapAverager.ROI.rowROIlastRow;
        colStart = handles.analysis.mapAverager.ROI.rowROIfirstCol;
        colEnd = handles.analysis.mapAverager.ROI.rowROIlastCol;
    case 'all'
        rowStart = handles.analysis.mapAverager.ROI.allROIfirstRow;
        rowEnd = handles.analysis.mapAverager.ROI.allROIlastRow;
        colStart = handles.analysis.mapAverager.ROI.allROIfirstCol;
        colEnd = handles.analysis.mapAverager.ROI.allROIlastCol;
    otherwise
end
roi = [rowStart rowEnd colStart colEnd];
% if strcmp('col', option)
%     Xvec = handles.analysis.mapAverager.xdata;
% elseif strcmp('row', option)
%     Xvec = handles.analysis.mapAverager.xdata;
% else
%     Xvec = [];
% end   
Xvec = [];


M = mapStack;
[r,c,p] = size(M);

% error catching for mapStack
if p < 2
    disp('mapROImeans: mapStack contains less than 2 maps; cannot evaluate means.');
    return
end
if sum(sum(sum(isinf(M))))>0
    M(isinf(M))=NaN;
end

% error catching for roi
if isempty(roi)
    roi = [1 size(M,1) 1 size(M,2)];
elseif numel(roi) ~=4
    disp('mapROImeans: ROI parameters incorrect.');
    return
end

% calculations
switch option
    case 'row'
%         for n = 1:p
%             m = M(roi(1):roi(2), roi(3):roi(4), n);
%             m = m';
%             V(n,:) = nanmean(m);
%         end
        for n = 1:p
            m = M(roi(1):roi(2), roi(3):roi(4), n);
            m = m';
            if roi(3) < roi(4)
                V(n,:) = nanmean(m);
            elseif roi(3) == roi(4)
                V(n,:) = m;
            end
        end
        vecX = roi(1):roi(2);
        vecY = V;
        vecAvg = nanmean(V);
        vecErr = nanstd(V)./sqrt(sum(~isnan(V)));
        titleTxt = 'average along rows';
        
    case 'col'
%         for n = 1:p
%             m = M(roi(1):roi(2), roi(3):roi(4), n);
%             %         m = m';
%             V(n,:) = nanmean(m);
%         end
        for n = 1:p
            m = M(roi(1):roi(2), roi(3):roi(4), n);
            %         m = m';
            if roi(1) < roi(2)
                V(n,:) = nanmean(m);
            elseif roi(1) == roi(2)
                V(n,:) = m;
            end
        end
        vecX = roi(3):roi(4);
        vecY = V;
        vecAvg = nanmean(V);
        vecErr = nanstd(V)./sqrt(sum(~isnan(V)));
        titleTxt = 'average along columns';
        
    case 'all'
        for n = 1:p
            m = M(roi(1):roi(2), roi(3):roi(4), n);
            V(n) = nanmean(nanmean(m));
        end
        vecX = 1;
        vecY = V;
        vecAvg = nanmean(V);
        vecErr = nanstd(V)./sqrt(sum(~isnan(V)));
        titleTxt = 'average over region';
end

if ~isempty(Xvec)
    if strcmp(option, 'row')
        vecX = -Xvec;
    elseif strcmp(option, 'col')
        vecX = Xvec;
    end
end

% put the output into the handles variable
handles.analysis.mapAverager.ROI.vecX = vecX;
handles.analysis.mapAverager.ROI.vecY = vecY;
handles.analysis.mapAverager.ROI.vecAvg = vecAvg;
handles.analysis.mapAverager.ROI.vecErr = vecErr;

% figure
if figFlag
    if get(handles.flipSignCheck, 'Value')
        S = -1;
    else
        S = 1;
    end
    
    hf = figure;
    set(hf, 'Units', 'normalized')
    pos = get(gcf, 'Position');
    set(hf, 'Position', [pos(1) pos(2) pos(3)*1.3 pos(4)*.7], 'Color', 'w');
    
    subplot(1,2,1);
    h = plot(vecX, S * vecY);
    pos = get(gca, 'Position');
    set(gca, 'Position', [pos(1) pos(2) pos(3) pos(4)*.85]);
    if strcmp(option, 'all')
        set(h, 'Marker', 'o');
    end
    grid on;
    axis tight;
    title('Individual expts');
    titleString = handles.analysis.mapAverager.currentGroup;
    ht = text(0, 0, titleString);
    set(ht, 'Units', 'normalized', 'Position', [-.2 1.2 0], ...
        'FontSize', 11, 'FontWeight', 'bold', ...
        'Interpreter', 'none');
    
    subplot(1,2,2);
    pos = get(gca, 'Position');
    set(gca, 'Position', [pos(1) pos(2) pos(3) pos(4)*.85]);
    h = errorbar(vecX, S * vecAvg, S * vecErr, 'bo-');
    axis tight;
    title('Avg & SEM');
    titleString = ['Analysis method: ' titleTxt];
    ht = text(0, 0, titleString);
    set(ht, 'Units', 'normalized', 'Position', [.2 1.2 0], ...
        'FontSize', 11, ...
        'Interpreter', 'none');
end

